

clear all
set more off
pause off
set matsize 11000


capture confirm file "Appendix_Results/"
if _rc mkdir "Appendix_Results/"


capture confirm file "Appendix_Results/gphs/"
if _rc mkdir "Appendix_Results/gphs/"



local region  Euro  

local response R_real



use "Data/Big_Pandemics.dta", clear
gen pandemic=1 
drop year
ren endyear year
sort year
merge 1:1 year using "Data/RealRateSince1310.dta"
sort year
replace pandemic=0 if pandemic==.
gen pandemic100k=cond(dead>=100000,pandemic,0)
cap drop Horizon
gen Horizon = _n-1 if _n<=41
tsset year


cap drop b0* se0* 
cap drop h0* l0* h1* l1*
tsset year
foreach x of local region {
	foreach y of local response {
	cap drop `x'`y'random
	cap drop `x'`y'c
	gen b0_`x'`y'=0 if Horizon==0
	gen se0_`x'`y'=0 if Horizon==0	
	gen h0_`x'`y'=0 if Horizon==0
	gen l0_`x'`y'=0 if Horizon==0
	gen h1_`x'`y'=0 if Horizon==0
	gen l1_`x'`y'=0 if Horizon==0

	
	cap drop raw
gen raw = .
replace raw = `x'`y'
/********** State space specification starts here *************/

constraint 1	[rstar]L.rstar		= 1		/* Random Walk */	
constraint 2	[raw]rstar			= 1		/* Definitional */


cap drop zero								/* generate zero line for figures */
gen zero = 0

tsset year
cap drop r_star

/************	 State space equations	 ******************/

sspace		(rstar L.rstar, state noconstant)		/// /* r-star */
			(raw rstar, noconstant ),				/// /* raw */
			covstate(identity) constraints(1/2) difficult /*method(kdiffuse)*/ iterate(25)
			
			qui{
			/* Generate smoothed estimates of r-star */
			predict r_star if e(sample), states smethod(smooth) equation(rstar)
			
			
	
			cap drop `x'`y'random
			gen	`x'`y'random = r_star

			}


								}
								
						}
						
						

*** get baseline LP

			
local region   Euro
foreach y of local response {
	foreach x of local region {
	eststo clear
		forv h=0/40 {
		
						cap drop  d`y'`h'
						gen d`y'`h' = f`h'.`x'`y'random- l.`x'`y'random
						reg d`y'`h' pandemic100k l(0/10).pandemic100k  l(2/11).`x'`y'random 
												
						eststo
						replace b0_`x'`y'  = _b[pandemic100k]                      if Horizon==`h'
						replace se0_`x'`y' = _se[pandemic100k]                     if Horizon==`h'
						replace h0_`x'`y'  = _b[pandemic100k]+2*_se[pandemic100k]  if Horizon==`h'
						replace l0_`x'`y'  = _b[pandemic100k]-2*_se[pandemic100k]  if Horizon==`h'
						replace h1_`x'`y'  = _b[pandemic100k]+1*_se[pandemic100k]  if Horizon==`h'
						replace l1_`x'`y'  = _b[pandemic100k]-1*_se[pandemic100k]  if Horizon==`h'
						}
						

		
	}
}

mkmat b0_EuroR_real , matrix(b) nomissing
mat bt = b'
matrix B_baseline = bt
		
						
						
*** placebo LP : random sets of pandemic events, 1000 replications

set seed 6789

qui forv r=1/1000 {
noi disp "placebo: replication `r' / 1000"

sum pandemic100k if pandemic100k==1
scalar Npandemics = r(N)
cap drop temp
gen temp = runiform() if pandemic100k ~= .
sort temp
cap drop placebo
gen placebo = cond(_n<=Npandemics,1,0)
sort year
replace pandemic100k=placebo

//noi sum placebo


local region  Euro    
local response R_real

foreach y of local response {
	foreach x of local region {
	eststo clear
		forv h=0/40 {
		
						cap drop  d`y'`h'
						gen d`y'`h' = f`h'.`x'`y'random- l.`x'`y'random
						reg d`y'`h' pandemic100k l(0/10).pandemic100k  l(2/11).`x'`y'random 
												
						eststo
						replace b0_`x'`y'  = _b[pandemic100k]                      if Horizon==`h'
						replace se0_`x'`y' = _se[pandemic100k]                     if Horizon==`h'
						replace h0_`x'`y'  = _b[pandemic100k]+2*_se[pandemic100k]  if Horizon==`h'
						replace l0_`x'`y'  = _b[pandemic100k]-2*_se[pandemic100k]  if Horizon==`h'
						replace h1_`x'`y'  = _b[pandemic100k]+1*_se[pandemic100k]  if Horizon==`h'
						replace l1_`x'`y'  = _b[pandemic100k]-1*_se[pandemic100k]  if Horizon==`h'
						}

	
}
}

mkmat b0_EuroR_real , matrix(b) nomissing
mat bt = b'
if `r'==1		matrix B = bt
	else 		matrix B = ( B \ bt )
				
}
	
mat list B

preserve

svmat B_baseline, names(b_baseline)

svmat B, names(b_)
svmat B, names(bsd_)
keep b* 
collapse (mean) b_* (sd) bsd_* 
gen obs=1
reshape long b_baseline b_ bsd_ , i(obs) j(h)
sum b*

gen u1 = b_+1*bsd_
gen l1 = b_-1*bsd_
gen u2 = b_+2*bsd_
gen l2 = b_-2*bsd_

gen Horizon = _n if _n<=41

		esttab , keep(pandemic*)
		cap drop Zero
		gen Zero=0
		twoway ///
		(  rarea u2 l2 Horizon if Horizon~=. , fc(blue)  fi(10) lcolor(white)  lpattern(none)  )  ///
		(  rarea u1 l1 Horizon if Horizon~=. , fc(blue)  fi(20) lcolor(white) lw(none)  lpattern(none)  )  ///
		(  line  Zero    Horizon if Horizon~=. , lc(gray)  lw(med) lp(dash)                     )  ///	
		(  line  b_     Horizon if Horizon~=. , lc(purple)  lw(thick) lp(dash)  )  ///
		(  line  b_baseline     Horizon if Horizon~=. , lc(blue)  lw(thick) lp(solid)  ///
				xlab(0(10)40, labsize(large)) ylab(, angle(0) labsize(large)  nogrid) xtitle("Years since pandemic event ended",  size(large)) ///
				ytitle("Percent",  size(large))  ///
				legend( on ring(0) pos(7) col(1) reg(col(none))  order(5 "Baseline" 4 "Placebo" ) ) )  ///
		, title()    ///
		xsize(5) ysize(4) scale(1.1) graphregion(color(white)) plotregion(color(white))  scheme(s1color)  ///
		name(`x'`y',replace)  ylab(-3(1)1)
		
graph save   Appendix_Results/gphs/FigureA7.gph, replace
graph export Appendix_Results/FigureA7.pdf, replace

restore
						
						
